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Abstract 

We study in detail the impact of dynamical quarks on the gluon mass generation mechanism, in 



the Landau gauge, for the case of a small number of quark families. As in earlier considerations, we 



assume that the main bulk of the unquenching corrections to the gluon propagator originates from 



£>. \ the fully dressed quark-loop diagram. The nonperturbative evaluation of this diagram provides 

the key relation that expresses the unquenched gluon propagator as a deviation from its quenched 

ON" 

l/"~) ■ counterpart. This relation is subsequently coupled to the integral equation that controls the mo- 

mentum evolution of the effective gluon mass, which contains a single adjustable parameter. The 

m ■ 

resulting nonlinear system is then treated numerically, yielding unique solutions for the modified 
gluon mass and the quenched gluon propagator, which fully confirm the picture put forth recently 

X- 

5_i , in several continuum and lattice studies. In particular, an infrared finite gluon propagator emerges, 

whose saturation point is considerably suppressed, due to a corresponding increase in the value 
of the gluon mass. This characteristic feature becomes more pronounced as the number of active 
quark families increases, and can be deduced from the infrared structure of the kernel entering in 
the gluon mass equation. 

PACS numbers: 12.38.Aw, 12.38.Lg, 14.70.Dj 



I. INTRODUCTION 



w 



The concept of a momentum- dependent gluon mass []J has been the focal point of con- 
siderable attention in recent years, mainly because it offers a natural and self-consistent 
explanation [2j for the infrared fmiteness of the (Landau gauge) gluon propagator and ghost 
dressing function, observed in large-volume lattice simulations, both in SU(2) ||3j and in 
577(3) [4|. Given the purely nonperturbative nature of the associated mass generation mech- 
anism, the Schwinger-Dyson equations (SDEs) constitute the most appropriate framework 
for studying such phenomenon in the continuum [5H8|. 

At the level of the SDEs, the general analysis finally boils down to the study of an 
integral equation, to be referred as the mass equation [9_|, |lOj , which governs the evolution 
of the dynamical gluon mass, m 2 (q 2 ), as a function of the momentum q 2 (see Eq. (\2.7\i 
below). It turns out that this special equation has been derived exactly, following the 
elaborate procedure explained in detail in 10] . This particular construction was carried 
out within the framework provided by the synthesis of the pinch technique (PT) [l|, llll-ll5J 



with the background field method (BFM) 16j, known in the literature as the PT-BFM 
scheme p, bj, |l8|. 

This exact mass equation, however, may not be treated in its complete form, because a 
particular field-theoretic ingredient entering in it is not fully known. Therefore, an approx- 
imate version of the full equation has been considered instead 10| , which depends on a free 
adjustable parameter. The detailed numerical study of this particular equation, for pure 
SU(3) Yang-Mills, revealed that its solutions depend strongly on the precise shape of the 
gluon propagators entering in its kernel. In particular, the q 2 region below 1 GeV 2 appears 
to be crucial for obtaining physically acceptable (i.e., positive-definite and monotonically 
decreasing) solutions. 



Interestingly enough, as has been shown in recent SDEs 19] and lattice studies 20], 
the gluon propagator, A(g 2 ), undergoes considerable changes in this particular region of 
momenta when a small number of light (dynamical) quarks is included in the theory, i.e., 
when "unquenching" takes place. This basic observation, in turn, motivates the further 
scrutiny of the mass equation under the special conditions that occur when the transition 
from pure Yang-Mills to real world QCD takes place. In particular, it is clearly of the 
outmost importance to establish, in quantitative detail, how the inclusion of dynamical 



quarks affects the entire mechanism of gluon mass generation (in the context of the so- 



called "scaling" solutions this issue was investigate in detail in 2l|, |22j). 

If one subscribes to the interpretation that the infrared finiteness of the lattice propagators 
is a consequence of the generation of such a mass, then the new unquenched lattice results 
of 20] would seem to indicate that the general mass generation picture persists in the 
presence of quark loops. In fact, the observed considerable suppression of the value of the 
saturation point of the unquenched propagators with respect to the quenched ones, would 
clearly suggest that the corresponding gluon mass increases [since A -1 (0) = m 2 (0)]. 

To explore this fundamental issue systematically, one needs a faithful description of the 
effects of the inclusion of quark-loops to the gluon propagator. To be sure, one could envisage 
the possibility of substituting the unquenched lattice data directly into the mass equation, 
in a spirit similar to that followed in the pure Yang-Mills case. However, we will refrain 
from using this latter set of data, and will resort instead to the full SDE approach presented 
in [19(, which allows the approximate inclusion of such effects in the quenched propagator. 

This particular approach assumes the presence of a small number of quark families, 
and treats the entire effect of "unquenching" as a "perturbation" to the quenched gluon 
propagator 1 . Then, the main bulk of the "unquenching" is attributed to the fully dressed 
quark-loop graph, while higher loop contributions are considered to be subleading. The 
inclusion of this particular nonperturbative loop induces considerable modifications to the 
shape of the resulting propagator, but does not provide precise information on the existence 
and value of the new saturation point. Indeed, such information may be only gathered in 
conjunction with the gluon mass equation, which was unavailable at the time of the original 



analysis. As a result, in the work of |19{| the existence of a saturation point was assumed, 
and its approximate value was obtained through simple extrapolation from the overall shape 
of the propagator in the region around 0.05 GeV 2 . 

Here, instead, our present knowledge of the mass equation allows us to carry out a 
thorough analysis of this particular dynamical question. Specifically, the equation describing 
the unquenching of the gluon propagator is coupled to to the mass equation, and the resulting 
non-linear system is solved numerically. In this way, one eventually obtains as a result the 



1 In the case of several quark families the changes induced may be rather profound, see, e.g., 23|, [2j], and 

references therein. 



unquenched gluon propagator in the entire range of physically relevant momenta, including 
the value of its saturation point. 

The main result of this article may be summarized as follows. The solution of the system 
reveals that the mass equation yields physical solutions in the presence of quarks, specifically 
for the typical cases Nj = 2 (two degenerate light quarks), and Nf = 2 + 1 + 1 (two degenerate 
light quarks and two heavy ones), considered in the recent lattice simulations [20|]. In fact, 
one establishes clearly that the value of the gluon mass at the origin increases as additional 
quarks are included in the theory. In addition, one obtains a concrete prediction for the 
quenched gluon propagator, which compares rather favorably with the aforementioned lattice 
results. We emphasize that throughout this analysis we make use of a single adjustable 
parameter, which is meant to model contributions to the mass equations stemming from the 
fully-dressed three gluon vertex. 

The article is organized as follows. In Section [Til we first define the basic quantities 
appearing in the problem, and introduce some important relations characteristic to the 
PT-BFM framework. Then, we describe the modifications that must be implemented at 
the level of the gluon SDE in order to consistently accommodate a dynamical gluon mass. 
Finally, we review the theoretical origin and properties of the two key dynamical equations, 
namely the gluon mass (integral) equation and the unquenching formula. Section IIHI is 
dedicated to the solution of the system formed when the aforementioned two equations are 
coupled to each other. To that end, we begin by reviewing the main ingredients appearing 
in them, and spell out the approximations employed in their treatment. Then we proceed to 
the numerical solution, describing in detail the iterative procedure adopted. The numerical 
results obtained are then compared with the latest lattice results, and some of the possible 
reasons for the observed deviations are discussed. Our conclusions are finally presented in 
Section [TVl 

II. GLUON MASS EQUATION WITH UNQUECHED GLUON PROPAGATORS 

The full gluon propagator (quenched or unquenched) in the Landau gauge assumes the 
general form 

iA^(q) = -zA(g 2 )P M! ,(g); P^(q) = 9^ - q^/q 2 , (2.1) 



where A(q 2 ) is related to the scalar form factor of the gluon self-energy U^q) = n(g 2 )P M;/ (g) 
through 

A-\q 2 ) = q 2 + tU(q 2 ). (2.2) 

In addition, the corresponding inverse gluon dressing function, J(q 2 ), is defined as 

A- V) = g 2 J(g 2 ), (2.3) 

and, as has been explained in [9|, |25|, it is intimately related with the definition of the 
renormalizat ion-group invariant QCD effective charge. 

As has been demonstrated in detail in the recent literature 



26 



271 ]. the nonperturbative 



dynamics of Yang-Mills theories gives rise to a dynamical (momentum-dependent) gluon 
mass, which accounts for the infrared finiteness of A(g 2 ). Specifically, from the kinematic 
point of view, we will describe the transition from a massless to a massive gluon propagator 
by carrying out the replacement A~ x (g 2 ) -» A~ x (g 2 ), where (Minkowski space) 

A~ m \q 2 ) = q 2 J m {q 2 )-m 2 {q 2 ). (2.4) 

Notice that the subscript "m" indicates that effectively one has now a mass inside the cor- 
responding expressions: for example, whereas perturbatively J{q 2 ) ~ lng 2 , after dynamical 
gluon mass generation has taken place, one has J m (q 2 ) ~ ln(g 2 + m 2 ). The presence of 
this mass, in turn, tames the Landau pole appearing in the (perturbative) effective charge, 



causing its saturation in the deep infrared 



28 



29|. 



From the dynamical point of view, the emergence of massive solutions from the gluon 
propagator SDE depends crucially on the inclusion of a set of special vertices, to be generi- 
cally denoted by V and called pole vertices [9] , which contain massless, longitudinally coupled 
poles. These vertices are added to the usual (fully dressed) vertices of the theory, and have a 
two- fold effect: they trigger the well-known Schwinger mechanism, thus endowing the gluon 
with a dynamical mass, while, at the same time, they guarantee that the Slavnov- Taylor 
identities (STIs) of the theory maintain exactly the same form before and after the mass 
has been dynamically generated. In particular, the transversality of Ii^ v {q) in the presence 
of masses is preserved, provided that one replaces all vertices Y appearing in the gluon SDE 
by 

r — > r = r m + v, (2.5) 



m 2 (q 2 ) 



1 1 
l + G{q 2 W 



g^x 






FIG. 1: The effective SDE satified by the dynamical gluon mass. The blue vertex corresponds to 
the vertex T m of Eq. (|2.5p , while the red vertex V indicates a pole vertex in which one of the gluon 
is leg is a background leg. Finally, the box singles out the quantity Y(k 2 ) which represents the 
purely two-loop dressed correction to the one-loop dressed mass equation kernel. 

with V having the precise structure that will make the new vertices T' satisfy the same 
formal STIs as the T before, but with the replacement A -1 — > A" 1 . 

It turns out that Jm{q 2 ) and m 2 (q 2 ) satisfy two separate, but coupled, integral equations 
of the generic type 2 [9| 

J m (q 2 ) = l+ JCi{k,q,m 2 ,A), 

Jk 



2/„2^ 



m (q 



JC 2 (k,q,m 2 , A), 



(2.6) 



such that /Ci, /C2 7^ 0, as q — > 0. In the equations above we have introduced the dimensional 
regularization measure J k = fi e Jd d k/(27r) d where d = 4 — e is the space-time dimension and 
H the 't Hooft mass. 

The exact closed form of the kernel JC2, and hence the full integral equation that deter- 
mines m 2 (q 2 ), has been derived in 10|, by appealing to the special properties of the pole 
vertices V. The resulting homogeneous equation is given by [l0|] (see also Fig. CE]) 

9 2 C A 1 



m (q ) 



1 + G(q 2 ) q 2 J k 



m 



2 {k 2 )^{k)A^{k + q)1C, u {k ) q) 1 



(2.7) 



with 



fC^(k t q) = [(k + q) 2 -k 2 ]{l-[Y(k + q) + Y(k)]}g, u 
+ [Y(k + q)-Y(k)}(q% u -2q tl q v ). 



(2.8) 



2 Note that in the rest of this article we will suppress the subscript "to" in A; it is understood that the 
gluon propagator is infrared finite. 



The quantity Y represents the subdiagram nested inside the two-loop dressed graph of Fig. [TJ 
given by 

Y(k 2 ) = ^Ak a jA a "(e)A^(e + k)r apP (-e-k,e,k), (2.9) 

with T ap s the full three-gluon vertex, and Ca the Casimir eigenvalue in the adjoint repre- 
sentation [C A = N for SU(N)]. 

Finally, the function G(q 2 ) corresponds to the g pu component of a special two point- 
function, which constitutes a crucial ingredient in a set of powerful identities, relating the 
conventional Green's functions to those of the BFM 30, |31j . For the case of the conventional 



gluon propagator, A, and the PT-BFM gluon propagator, denoted by A, the relevant identity 
reads 

A(g 2 ) = [l + G(g 2 )] 2 A(g 2 ); (2.10) 

its application at the level of the corresponding SDEs leads to the appearance of the factor 
1 + G(q 2 ) inEq. (I2T7D. 

The numerical treatment of Eq. (12. 7p . under certain simplifying assumptions regarding 
the form of Y(k 2 ) (see next section), gave rise to solutions that display the basic qualita- 
tive features expected on general field-theoretic considerations, and employed in numerous 
phenomenological studies. In particular, the m 2 (q 2 ) obtained are monotonically decreasing 
functions of the momentum and vanish rather rapidly in the ultraviolet [l|, Il0 |. 

As explained in the Introduction, the purpose of the present work is to study the effects 
induced on the gluon mass by the inclusion of a small number of light quark families. To 
that end, one needs to solve Eq. (12.7ft using unquenched gluon propagators, A N (q 2 ), namely 
carrying out the substitution A(g 2 ) — > A N Aq 2 ) inside the integral of the r.h.s.. 

The corresponding expression for A N (q 2 ) will be obtained following the SDE approach 



presented in [19j. Specifically, one assumes that the main bulk of the unquenching effect is 
captured by the (fully dressed) one-loop diagram of Fig. [2] (b), neglecting, at this level of 
approximation, all contributions stemming from (higher order) diagrams containing nested 
quark loops. 

Denoting the contribution of this special diagram by 

X^(q)=X(q 2 )P^(q), (2.11) 

we have that 

X(q 2 ) = - 9 j I Tr [-fSihfyk, -k - q, q)S(k + q)] . (2.12) 



A" 



^ 2 ) = [1 + G(q 2 )]- 2 X(q 2 ) + m 2 Nf (q 2 )-m 2 (q 2 ) + A_V) 
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FIG. 2: Schematic representation of the unquenched propagator (a), corresponding to Eq. (|2.15p . 
and some of the ingredients [(b), (c), and (d)] entering in it. In particular, (6) represents the quark 
loop, which, in the approximation employed, is the only source of quark-dependence. The quark 
propagators entering in (6) are solutions of the gap equation depicted in (c), where too denotes the 
appropriate current mass. Finally, the function G(q 2 ) is obtained from the relation shown in (d), 
namely Eq. f|3.2j> . The quantities obtained from the lattice are also indicated. Note that the term 
X 2 (q 2 ) will be determined dynamically, once the mass equation of Fig. [T] is coupled to (a). 



where S denotes the full quark propagator, and the vertex T^ corresponds to the PT-BFM 
quark-gluon vertex, satisfying the QED-like WI 



iq'T^k, -k -q,q) = S~\k) - S~\k + q). 



(2.13) 



Of course, in the case of including various quark loops, corresponding to different quark 
flavors Nf, the term X^ v (q) in Eq. ( 12. lip is replaced simply by the sum over all quark loops, 



i.e., 



X^(g)^]Txf(g). 
/ 



(2.14) 



Then, the detailed analysis of |l9| reveals that the unquenched gluon propagator A N (q 2 ) 
may be expressed as a deviation from the quenched propagator A(q 2 ), namely (Euclidean 



space) 



A(q 2 



A^(g 2 ) = ^ ^- , , (2.15) 

1 + |X(g 2 ) [1 + G(g 2 )]- 2 + A 2 (g 2 ) j A(g 2 ) 

where the quantity 

\ 2 (q 2 )=m 2 Nf (q 2 )-m 2 (q 2 ), (2.16) 

measures the difference induced to the gluon mass due to the inclusion of quarks. 



An important point, which can be established formally and confirmed numerically 19], 
is that the nonperturbative X(q 2 ) vanishes at the origin, X(0) = 0, exactly as it happens 
in perturbation theory. This implies that the mass equation is not directly affected by the 
presence of dynamical quarks, as the methodology used to derive Eqs. (j2.7p . ( 12. 8 p and (12. 9p 
is left invariant by the unquenching procedure. However, the solutions of the corresponding 
equation will be different from those obtained in the quenched case, as the kernel /C MI/ will 
change, due to the modifications induced by X(q 2 ) in the overall shape of the propagator 
throughout the entire range of momenta. Thus, the inclusion of quark loops affects the 
value of the saturation point of the gluon propagator, not directly through the presence of 
the X(q 2 ), rather indirectly through the generation of a non vanishing mass difference A 2 (g 2 ). 

The next nontrivial step is therefore to treat the mass equation Eq. ( j2 . Tf) and the un- 
quenching master formula of Eq. (12. 15ft as a coupled system, and determine both m 2 N and 
Ajv Aq 2 ). This rather involved task will be undertaken in detail in the next section. 

III. NUMERICAL ANALYSIS 

Before entering into the technical details of the solution of the system, it is important 
to comment on a fundamental qualitative difference between the present situation and the 



treatment followed in 10j. Specifically, as already explained in the previous section, even in 
the absence of quarks, Eq. (12.71) forms part of a system of coupled equations, namely that of 
Eq. ( 12. 6p . However, given that the corresponding equation for J m (q 2 ) is only partially known, 
the approach adopted in [l0( was to treat the gluon propagators appearing in Eq. ( 12. 6ft as 
external quantities, using the lattice results of [4| for their functional form. As a consequence, 
the intrinsically non-linear equation (I2.6P (recall Eq. ( I2.4p ) was converted into a linear one. 
As a result of this linearity, one obtained a continuous family of solutions, parametrized 
by the values of a multiplicative constant; the actual solution chosen was the one that 



reproduced the saturation point of the lattice gluon propagator that was used as input in 
Eq. d23J. 

Here, instead, even though we still decouple the equation for J m (q 2 ), the mass equation 
retains its non-linear nature due to the form of the unquenched gluon propagator, A N (q 2 ), 
that will be inserted in it. Specifically, as can be appreciated from Eq. ( I2.15p . the X 2 (q 2 ) 
appearing in the denominator of A N (q 2 ) introduces a non-linear term into Eq. ( 12. 6p . Con- 
sequently, the solution obtained will be unique, and the saturation point of the resulting 
A N (q 2 ) will emerge as a prediction of the theory rather than as an input from the lattice. 

A. Main ingredients and basic assumptions 

We next proceed to a brief description of the ingredients entering into the two equa- 
tions (12. 6p and (I2.15P composing the system under study, together with an account of the 
most important underlying assumptions. 

(i) The determination of the quantity X(q 2 ) is of central importance, since it approx- 
imates the entire effect of the inclusion of the active quarks. In order to evaluate it from 
Eq. (I2.12p . one needs specific expressions for the full quark propagator S and the fully- 
dressed vertex I\ Evidently, the total X(q 2 ) emerges as the sum of the individual flavor 
contributions, in accordance with Eq. (j2.14p . 

The nonperturbative form of each quark propagator entering into X(q 2 ) is obtained from 



the standard gap equation [32] , supplemented by an appropriate current mass term, in order 



to make contact with the lattice results of 20j. In this latter simulation, the gluon (and 



ghost) propagators have been evaluated from large volume configurations (up to 3 3 x 6 
[fm 4 ]), generated from a lattice action that included (twisted mass) fermions. Specifically, 
one employed two light degenerate quarks (Nf = 2), with a current mass ranging from 20 to 
50 MeV, or two light and two heavy quarks (Nf = 2 + 1 + 1), with a strange (charm) quark 
current mass roughly set to 95 MeV (1.51 GeV). 

For the quark propagator, S(p), we use the standard parametrization 

S-^p) = -i \A( v )i - B(p)\ = -iA(p) \i - M(p)\ , (3.1) 

where A(p) represents the "quark wave-function", while the ratio Ai(p) = B(p)/A(p) de- 
notes the dynamical quark mass, and extracts them from the corresponding system of inte- 
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FIG. 3: (color online) The inverse of the quark wave- functions (top), and dynamical quark masses 
(bottom), obtained from the quark gap equation for three different values of the current mass: 
mo = 41.2 MeV (black, continuous), too = 95 MeV (red, dotted) and too = 1.51 GeV (blue, 
dashed). 

gral equations. The inverse of the quark wave-functions, 1/A(p), and the dynamical quark 
masses, Ai(p), used in the present work, are shown in the top and bottom panels of Fig. |3l 
respectively. They have been obtained from the gap equation of [32[, where the quenched 
lattice data of |4j were used as input. 

As can be seen from Fig. [3] the two degenerate light quarks acquire a physical mass around 
300 MeV, while the two heavier ones get physical masses around 440 MeV (strange) and 2.2 
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FIG. 4: The full nonperturbative quark loop contribution X(q 2 ) for the two cases Nf = 2 (black, 
continuous) and Nf = 2 + 1 + 1 (blue, dashed). 



c masses are consistent 
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GeV (charm). It should be stressed that these values for the quar 
with the ones generally employed in phenomenological calculations 

As for the PT-BFM quark-gluon vertex T, the fact that it satisfies the Abelian WI of 
Eq. ( I2.13P allows one to model it by means of the standard Abelian Ansatze existing in 



the literature 



35 



361 ] . In particular, T is expressed entirely in terms of the quantities A(p) 



and A4(p), with no reference to auxiliary ghost Green's functions, a fact that simplifies 
considerably the task of computing X(q 2 ). 

When all the aforementioned ingredients are put together, the two different X(q 2 ), cor- 
responding to the cases Nf = 2 (black, continuous) and Nf = 2 + 1 + 1 (blue, dashed), are 
shown in Fig. |U 



Note finally that, as one can establish formally and confirm numerically 19), the non- 
perturbative X(q 2 ) vanishes at the origin, X(0) = 0, exactly as it happens in perturbation 
theory. Evidently, the inclusion of quark loops does not affect directly the value of the sat- 
uration point of the gluon propagator. However, as already explained, the saturation point 
will be indirectly affected, due to the modifications induced by X(q 2 ) in the overall shape of 
the propagator, throughout the entire range of momenta. 

(ii) As already mentioned, the quenched lattice data for A(q 2 ) [4| constitute the starting 
point for obtaining from Eq. (J2.15P a prediction for A N (q 2 ). In addition, in the Landau 
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0.0 0.2 0.4 0.6 0.8 1.0 



FIG. 5: The curves formed by the set of the pairs (Ca s , a s ) yielding physical solutions to the full 
mass equation (|2.7p . for the various numbers of quark families Nf. quenched (red, dots), Nf = 2 
(green, triangles) and Nf = 2 + 1 + 1 (blue, squared). At the values of C indicated by the stars, 
the solutions obtained corresponds to the expected values of the strong coupling a s at the scale 
chosen (p = 4.3 GeV); one has C = 9.04 for N f = (a s = 0.22), C = 8.48 for N f = 2 (a s = 0.285), 
and C = 8.64 for N* = 2 + 1 + 1 (a s = 0.33). Finally, the values of a s obtained for the case in 
which Y is simply kept at the lowest order perturbative value (i.e., C = 1) are a s = 0.93 (Nf = 0), 
a s = 1.09 (N f = 2), and a s = 1.25 (N f = 2 + 1 + 1). 

gauge, the quantity 1 + G(q 2 ), appearing in both Eos. (12.61) and f)2.15p . is linked to the 
inverse of the ghost dressing function F(q 2 ) through 37], 



F- L (q 2 )^l + G(q 2 ). 



(3.2) 



This relation, which is valid to a very good approximation, and becomes an exact equality 
at q 2 = 0, allows one to use the lattice results of j4j for the ghost dressing function, in order 
to determine G(q 2 ). Notice that both sets of data will be renormalized at the scale /i = 4.3 
GeV, which is the last available point in the ultraviolet. At this scale, the expected value of 
the (quenched) strong effective charge is a s = g 2 /4:ir = 0.22 |39l. |40|. 

(in) The function Y(k 2 ) represents a crucial ingredient of Eq. (12. 7\i . However, its exact 
closed form is not available, mainly because our present knowledge of the full three-gluon 
vertex, entering in its definition, is relatively limited; we must therefore resort to suitable 
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approximations. In particular, we will employ the lowest-order perturbative expression 
for Y(k 2 ), obtained from Eq. (12. 9p by substituting the tree-level values for all quantities 
appearing there. Within this approximation, and after carrying out momentum subtraction 



renormalization (MOM) at k 2 = /i 2 , one obtains 10] 



^ 2 > = -^ 1o 4 < 3 - 3 > 

where a s is the value of the coupling at the subtraction point chosen. This simple approx- 
imation will be compensated, in part, by multiplying Yr(/c 2 ) by and arbitrary constant C, 
i.e., by implementing the replacement Yr(A; 2 ) — > C 1r(/c 2 ) treating C as a free parameter. In 
this heuristic way, one hopes to model further corrections that may be added to the "skele- 
ton" result provided by Eq. (13.31) . This particular procedure has been proved sufficient for 
obtaining in the quenched case physically meaningful solutions for m 2 (q 2 ) for reasonable val- 



ues of the constant C lOj. As the numerical study presented in Fig. shows, this continues 

to be true also in the unquenched case. 

(iv) A particularly helpful constraint may be obtained by taking the q 2 — > limit of the 

mass equation (12.71) . which gives (passing to spherical coordinates, and setting y = k 2 ) 

3 f°° < -\' 

m 2( ) = -— as C A F(0)J dym 2 (y))C Nf (y)- K Nf {y) = [[1 -2CY(y)]y 2 A 2 Nf (y)} , 

(3.4) 

where the prime denotes derivatives with respect to y. The usefulness of this relation lies in 

the fact that, when used in combination with the unquenched propagators obtained following 

the approach of [19J], it allows one to deduce how the running mass will be affected by the 

presence of dynamical fermions. 

To see how this works in detail, consider first the quenched case; the kernel /Co, con- 
structed from the corresponding lattice propagator data, is shown on the top panel of Fig. El 
As can be appreciated there, for increasing values of C one observes the appearance of a 
negative area, in the region of momenta q 2 < 1 GeV 2 . It is precisely the contributions 
from this negative region that counteract the minus sign in front of the mass equation ( 12. 7p . 
thus making possible the existence of positive-definite and monotonically decreasing solu- 
tions (l0|. 

Let us now repeat this exercise for the unquenched case with Nf = 2+1+1. On the bottom 
panel of Fig. [6] we show the kernel /C 2 +n-ij the A 2+ i+i(g 2 ) used to construct it is obtained 
from Eq. (J2.15J) . using the X that corresponds to the aforementioned quark configuration, 
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FIG. 6: (color online) Comparison between the kernel /Cjy.(y) of Eq. (|3.4p for the quenched case 
(top), and the unquenched case with Nf = 2 + 1 + 1 (bottom). The plots are drawn with the same 
scale to facilitate the comparison. 

and setting A 2 = 0. As we will see in the next section, this particular propagator emerges 
precisely after the first step of the iterative procedure employed to solve the SDE system. 
Note also that, due to the decoupling of the heavy fermions [19|, the kernel /C 2 , corresponding 
to the case Nf = 2, is practically identical to /C2+1+1. 

As it is evident from Fig. [61 the addition of dynamical fermions affects the shape of 
the kernel significantly, as it now displays a much shallower negative region compared to 
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the quenched case. Once this kernel is inserted into the mass equation (for the successive 
iterative steps), one then observes that in order to keep satisfying the constraint Eq. (13. 4p 
with a positive m 2 (0), the dynamical mass must increase (and therefore the saturation point 
of the unquenched propagator decrease), in order to compensate for the suppressed negative 
region. This increase is counteracted in part by the fact that the coupling constant a s also 
rises when adding quarks [at [i = 4.3 GeV, one has a 30% increase for the Nf = 2 case (from 
0.22 to 0.285), and an additional 20% increase (from 0.285 to 0.33) for the N f = 2 + 1 + 1 
easel . 

n 

This qualitative conclusion is supported by both the SDE study of [19] as well as the 



lattice results of 20j, and will be fully confirmed through the explicit numerical computations 
carried out in the next subsection (see in particular Fig. [TD], where the final running masses 
for the different cases are presented). 

B. Solving the system and comparing with the lattice 

The system of SDEs that we consider is composed of Eqs. (12. 7\i . (12.121) and (I2.15p . supple- 
mented by the quark gap equation. The initial condition is provided by the quenched SU(3) 
gluon propagator and ghost dressing function obtained in the lattice simulations of [4j , which 
will be also used to determine the initial values of the form factors A(p) and Ai(p). All 
calculations will be performed at fi = 4.3. 

The algorithm that we adopt consists of the following main steps. 

(z) Using the quenched propagator as an input of the first iterative step, one determines 
the quark form factors A{j>) and .M(p) by solving the quark gap equation; (ii) A{j>) and 
Ai(p) are substituted into Eq. (12.121) . and the corresponding value of the quark loop diagram 
X(q 2 ) is evaluated; (Hi) The preliminary form of A N (q 2 ) is determined from (12. 15ft . employ- 
ing initially A 2 (g 2 ) = 0, with the quenched mass m 2 (q 2 ) obtained from the solution of the 
mass equation (12.71) corresponding to the quenched lattice propagator; (iv) The unquenched 
propagator A N Aq 2 ) of the previous step is substituted into the mass equation (12. 7p in order 
to determine the associated unquenched dynamical gluon mass m 2 N (q 2 ), and therefore the 
corresponding A 2 (g 2 ); (v) At this point the latter quantity is inserted back into the master 
equation (12.151) . and the loop starts again, until convergence, determined by the stability of 
the quantities involved, has been reached. 
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FIG. 7: (color online) Convergence of the iterative procedure for solving the SDE system in the 
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2 + 1 + 1 case. Top: unquenched propagator A N (q 2 ); Bottom: mass difference X 2 (q z 



In Fig. [7] we show how the procedure described above converges rather rapidly to a 
stable solution for both A N (q 2 ) (top panel) as well as A 2 (g 2 ) (bottom panel). In particular, 
observe that at the first iterative step one has A 2 (g 2 ) = 0; even though this gives rise to 
a propagator with the same saturation point as in the quenched case (red circles curve on 
the top panel), the presence of the fermion loops alter in a rather marked way the overall 
shape of the function. As discussed in the previous section, this affects in turn the kernel of 
the mass equation, resulting in a A 2 (g 2 ) > already at the second step, which then forces 
the unquenched propagator to have a lower saturation point with respect to the quenched 
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FIG. 8: (color online) The unquenced gluon propagators obtained from our analysis, for Nf = 2 
(upper panel) and Nf = 2 + I + I (lower panel), compared with the lattice data of [20] for the same 



cases. 



case (green stars curves on both panels). Notice that, as anticipated, A 2 (g 2 ) increases at 
each successive iteration, resulting in an unquenched propagator that is suppressed in the 
infrared with respect to the quenched case. 

In Fig. [8] we present the central result of this analysis. In particular, we plot the prop- 
agators obtained when convergence of the above mentioned iteration procedure has been 
reached, and compare them with the corresponding unquenched lattice data, recently re- 
ported in [20| . We observe a rather good agreement between our theoretical predictions and 
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FIG. 9: (color online) The gluon dressing functions, q 2 A N (q 2 ), corresponding to the propagators 
of the previous figure. 

the lattice computation for both values of Nf, for the available range of physical momenta. 
A notable exception to this fair coincidence between curves is the saturation point of the 
Nf = 2 + 1 + 1 case; specifically, the value obtained from our SDE analysis is 20% higher 
than that found in lattice simulations. Similar conclusions can be drawn by observing the 
plot corresponding to the gluon dressing function q 2 A N (q 2 ) shown in Fig. ED while one has 
an excellent agreement in the case of two degenerate light quarks, when two heavier quarks 
are added the SDE solution tends to mildly overestimate the amplitude of the characteristic 
peak located in the intermediate momentum region. 

The corresponding dynamical gluon masses, m\(q 2 ) and m 2 +1+1 (q 2 ), are shown in Fig.[T0~| 
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FIG. 10: (color online) Solution of the mass equation yielding the dynamically generetaed gluon 
mass for Nf = 2 (red dotted line) and Nf = 2 + 1 + 1 (blue dashed line). In the deep infrared one 
has m2(0) = 413 MeV, and 7772+1+1(0) = 425 MeV. For comparison we also show the quenched 
gluon mass (black continuous line) obtained from the quenched lattice propagator, in which case 
m (0) = 376 MeV. 

for comparison, we also plot the quenched solution, m 2 (q 2 ), obtained from Eq. (12. 6p when 
the quenched lattice propagators of [4] are used as input. In particular, the corresponding 
saturation points give 7712(0) = 413 MeV and 7712+1+1(0) = 425 MeV (at \x = 4.3 GeV), which 
should be compared with the value m(0) = 376 MeV obtained for the quenched case. The 
results captured in Fig. [10] are particularly important, because they demonstrate clearly 
that the mass generation mechanism established for pure Yang-Mills continues to operate 
in QCD-like circumstances. 

Let us finally comment on some possible factors that might account for the aforementioned 
discrepancies between our predictions and the lattice data for Nf = 2 + 1 + 1. 

On the one hand, it is generally known that the deep infrared region of Green's functions 
simulated on the lattice is affected by volume artifacts. As a result, one may expect moderate 
fluctuations in the value of A w (0), as the lattice volume (spacing) is increased (decreased). 
Notice also that the access to large volumes and small lattice spacings is restricted by the 
currently available unquenched gauge configurations; for example, in the lattice simulations 
of [20|], one has, as already said, a volume corresponding to at most 3 3 x 6 fm 4 . On the other 



20 



hand, it is reasonable to assume that the observed deviation might signal the mild violation 
of one of the assumptions underlying the derivation of the unquenching formula (12.151) . In 
particular, it is natural to expect that one of our basic operating assumption, namely that 
the quark-loop contributions constitute a "perturbation" of the quenched propagator, will 
become progressively less accurate as the number of active flavors increases. It is therefore 
possible that from Nf > 2 onward we begin to perceive the onset of additional effects, 
not captured by (j2.15p . For instance, the identification of the "pure" gluonic diagrams of 
the gluon SDE with the gluon propagator A(g 2 ) obtained in quenched lattice simulations, 
may no longer be entirely reliable, due to the presence of (increasingly appreciable) effects, 
coming from higher-order quark sub diagrams. 

Clearly, an additional theoretical uncertainty originates from the approximate (pertur- 
bative) treatment of the quantity Y(k 2 ). In particular, the parameter C may only model, 
to some extent, unknown contributions that display a logarithm momentum dependence, 
as in Eq. ( 13.3}) . but cannot account for terms with a different functional form. Moreover, 
the use of an Ansatz for the vertex T^ entering into the definition of X(q 2 ) may induce 
further error, due to the fact that its transverse (automatically conserved part) is in general 
undetermined. 

IV. CONCLUSIONS 

In this work we addressed the basic theoretical question of whether the mass generation 
mechanism established in pure Yang-Mills studies persists in real-world QCD. Specifically, 
by employing a methodology relying mainly on the SDEs that describe the gluon two-point 
sector within the PT-BFM framework, we studied in quantitative detail how the inclusion 
of dynamical quarks affects the generation of the momentum-dependent gluon mass, in the 
Landau gauge. This important issue is especially relevant and timely, given the qualitative 
picture that appears to emerge from recent unquenched lattice simulations. Our results 
demonstrate clearly that the gluon propagator still saturates in the infrared (therefore, a 
gluon mass is indeed generated), and that the saturation point is progressively suppressed, 
as the number of quark flavors increases. 

Several of the salient dynamical features pertinent to the transition from the quenched 
to the unquenched theory have been presented in the preliminary study of [19]. The main 
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novelty of the present analysis resides in the fact that, unlike [19J . where the saturation point 
of the gluon propagator was estimated by means of an extrapolation procedure, here it is 
determined explicitly from the solution of the gluon mass equation. Therefore, the results 
obtained constitute a genuine prediction of the formalism employed, emerging entirely from 
the inherent dynamical equations. 

It is worth pointing out that the ingredients obtained from the present analysis may 
be used to determine the flavor-dependence of the QCD effective charge, a s (q 2 ). Specif- 
ically, as has been shown in 25], this latter quantity may be defined in terms of 
the functions J(q 2 ) and G(q 2 ), through the renormalization-group invariant combination 
Q-s 1 ^ 2 ) = a 7 1 (/ i2 )[l + G(q 2 ; ^ 2 )} 2 Jm(<f'i (i 2 ) where /i is the renormalization (subtraction) 
point chosen, within an appropriate renormalization scheme. The direct determination of 
J(q 2 ) from its own dynamical equation, shown schematically in Eq. 02 .6j) , is thwarted from 
the fact that the main ingredients entering in it are the fully dressed three-gluon vertex BQ 2 



(studied in [4l|]) and the (largely unexplored) four-gluon vertex BQ 3 , for arbitrary values 
of their momenta. Instead, one may employ Eq. ( 12 Ah to compute J{q 2 ) indirectly, using a 
combined approach based on the knowledge of the full propagator from the lattice, and the 
corresponding gluon mass from solving Eq. (12 .7p . (For the possibility of the direct extraction 
of the quantity J m (q 2 ; jjl 2 ) = [1 + G(q 2 ; /i 2 )] 2 J m (g 2 ; /i 2 ) from specialized lattice simulations, 



see 42J). Work in this direction is already in progress. 

Finally, it would clearly be important to study the infrared dynamics of the gluon prop- 
agator, quenched and unquenched, away from the Landau gauge. The SDEs formulated 
within the PT-BFM framework provide a natural starting point for such an investigation in 
the continuum. In fact, it would be particularly relevant to determine whether the detailed 
mass generation mechanism established in the Landau gauge persists for other values of the 
gauge-fixing parameter, such as, for example, the Feynman gauge, which is known to be of 
central importance for the PT. To be sure, further reliable input from the lattice would be 
invaluable for accomplishing such a demanding task, given that, even though the formal pro- 
cedure for implementing covariant R^ gauges on the lattice has been already reported 



43), 



the available simulations are restricted to modest size volumes, and ^ < 1 only AA\. We 
hope to report progress in some of these directions in the near future. 
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